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ABSTRACT 

The propagation of high-energy cosmic rays through giant molecular clouds constitutes a 
fundamental process in astronomy and astrophysics. The diffusion of cosmic-rays through these 
magnetically turbulent environments is often studied through the use of energy-dependent diffu- 
sion coefficients, although these are not always well motivated theoretically. Now, however, it is 
feasible to perform detailed numerical simulations of the diffusion process computationally. While 
the general problem depends upon both the field structure and particle energy, the analysis may 
be greatly simplified by dimensionless analysis. That is, for a specified purely turbulent field, the 
analysis depends almost exclusively on a single parameter - the ratio of the maximum wavelength 
of the turbulent field cells to the particle gyration radius. For turbulent magnetic fluctuations 
superimposed over an underlying uniform magnetic field, particle diffusion depends on a second 
dimensionless parameter that characterizes the ratio of the turbulent to uniform magnetic field 
energy densities. We consider both of these possibilities and parametrize our results to provide 
simple quantitative expressions that suitably characterize the diffusion process within molecular 
cloud environments. Doing so, we find that the simple scaling laws often invoked by the high- 
energy astrophysics community to model cosmic-ray diffusion through such regions appear to be 
fairly robust for the case of a uniform magnetic field with a strong turbulent component, but 
are only valid up to ~ 50 TeV particle energies for a purely turbulent field. These results have 
important consequences for the analysis of cosmic-ray processes based on TeV emission spectra 
associated with dense molecular clouds. 

Subject headings: Cosmic Rays - diffusion - ISM - molecular clouds 
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1. Introduction 

Observations of 7-rays associated with regions 
of dense molecular gas provide important clues 
about how cosmic-rays (CR's) are injected within 
our galaxy. However, a proper treatment of this 
problem requires an understanding of how CR's 
diffuse through turbulent environments. While 
this subject has received considerable attention 
since the pioneering works of Jokipii (1966) and 
Kulsrud & Pearce (1969), the exact nature of par- 
ticle transport remains unresolved. 

A standard approach to the problem invokes 
the use of the spherically symmetric diffusion 
equation 



dt 



D_d_ 

R 2 OR 
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dR + dE r 



(Pf) + Q, (1) 



where / = f(E p , R, t) is the distribution of parti- 
cles as a function of energy, distance, and time; 
P = —(dEp/dt) is the continuous energy loss 
rate; Q = Q(E p ,R,t) is the source function; and 
D = D(Ep) is the energy-dependent diffusion coef- 
ficient. A simplified solution to this diffusion equa- 
tion may be obtained by assuming a power-law 
injection spectrum, /; n j oc E~ a , and a power-law 
diffusion coefficient, 



D(E P ) = D 
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E n 



lOGeV 



(2) 



in the energy regime where t pp is independent 
of energy (we note that values of S = 1/2 and 
Dio <~ 10 26 ~ 28 cm 2 s _1 are typically assumed for 
molecular cloud environments — see, e.g., Aharo- 
nian & Atoyan 1996; Torres et al. 2003; Gabici 
et al. 2009). As shown by Aharonian and Atoyan 
(1996), the solution to the diffusion equation in 
such a case can be approximated as: 



f(E p ,R,t) 
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-Rdiff = Rdiff (E p ,t) — 
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is the "diffusion radius" corresponding to the ra- 
dius of the sphere out to which particles with en- 
ergy E p effectively propagate after a time t. In the 
limit that t <C t pp , the "diffusion radius" simplifies 
to R dis = 2y/D(E p )t. 

In this paper, we investigate how high-energy 
CR's propagate through molecular cloud-like envi- 
ronments by instead using a modified numerically 
based formalism developed for the general study 
of cosmic-ray diffusion by Giacalone & Jokipii 
(1994). This formalism has already been used to 
study the transport of cosmic rays in chaotic mag- 
netic fields with Kolmogorov turbulence (Casse et 
al. 2002) and has been applied successfully in sev- 
eral specific contexts (see, e.g., Kowalenko & Melia 
1999; Casse et al. 2002; De Marco et al. 2007; 
Wommer et al. 2008; Fraschetti & Melia 2008). 

The first goal of this work is to extend the gen- 
eral treatment of Casse et al. (2002) by explor- 
ing a greater dynamic range of wavelengths over 
which turbulence acts and by considering Kraich- 
nan, Bohm and Kolomogorov turbulence for two 
magnetic field configurations: 1) a purely turbu- 
lent field; and 2) a uniform magnetic field with 
a strong turbulent component. The second goal 
of this work is to provide a baseline analysis for 
the propagation of <~ 1 — 10 TeV cosmic-rays in 
molecular cloud environments. 

As we shall sec, CR diffusion in purely turbulent 
fields depends primarily on a single dimensionless 
parameter 

T _ Amax 

A max — „ , \o) 

Rg 

where A max represents the longest turbulent field 
wavelength and R g is the particle gyration radius 
in a uniform field of the same magnetic energy 
density as that of the turbulent field. This param- 
eter is related to the particle rigidity p through 
the expression A max = 2-k/ p. In the second case, 
CR diffusion also depends on a second dimension- 
less parameter — the ratio of turbulent field energy 
density to the uniform field energy density. As we 
shall see, the result of our work indicates that the 
diffusion coefficients often invoked to describe CR 
diffusion through molecular cloud environments 
appear to be valid for < 50 TeV cosmic rays prop- 
agating in a purely turbulent field, and appear to 
be fairly robust for the case of a uniform magnetic 
field with a strong turbulent component. 
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Our paper is organized as follows. The relevant 
properties of molecular clouds are briefly reviewed 
in §2, where we also outline our treatment of these 
environments. The scheme for generating the tur- 
bulent magnetic field is presented in §3, and the 
equations that govern the motion of CR's are di- 
mensionalizcd in §4. Solutions to these equations 
are presented in §5 for purely turbulent fields, and 
in §6 for a uniform field with a strong turbulent 
component. We compare and contrast the results 
of our work to those of Casse et al. (2002) in 
§7. We then consider what effects our results have 
on previous treatments of CR diffusion through 
molecular clouds in §8, and summarize our work 
in §9. 

2. Giant Molecular Cloud Environments 

Typical giant molecular clouds (GMCs) con- 
tain a total mass of ~10 5 M Q within physical size 
scales of tens of parsecs, and, as such, have mean 
densities of n# 2 ^100 cm~ 3 . However, these large 
complexes are highly nonuniform, exhibiting hi- 
erarchical structure that can be characterized in 
terms of clumps (R~l pc, n# 2 ^10 3 cm~ 3 ) and 
dense cores (i?~0.1 pc, n# 2 ^10 4 -10 5 cm~ 3 ) sur- 
rounded by an interclump gas of density hh 2 ~ 5- 
25 cm~ 3 . 

Exactly how the magnetic field is partitioned 
within GMCs is not yet known. In the simplest 
case, where flux freezing applies, the magnetic field 
strength B in the interstellar medium would scale 

1/2 

with the gas density nn 2 according toBot . It 
is noteworthy, then, that an analysis of magnetic 
field strengths measured in molecular clouds yields 
a relation between B and nn 2 of the form 

b - 10 " g (to^)"' 4 '' < 6 » 

though with a significant amount of scatter in the 
data used to produce this fit (Crutcher 1999; but 
see also Basu 2000). This result is consistent with 
the idea that nonthermal lincwidths, measured to 
be ~1 km s _1 throughout the cloud environment 
(e.g., Lada et al. 1991), arise from MHD fluctua- 
tions. 

The exact nature of the magnetic turbulence 
is not well-constrained, although magnetic fluc- 
tuations are typically assumed to have a power- 
law spectrum such that their intensity at a given 



wavenumber scales according to (5Bk) 2 ~ fc~ r , 
with indices typically taken to be T = 1 (Bohm), 
T = 3/2 (Kraichnan) or T = 5/3 (Kolmogorov). 
In addition, the range in wavelengths over which 
these fluctuations occurs is not well known, al- 
though it is reasonable to assume that the up- 
per end corresponds to the lengthscale over which 
the fluctuations are generated. (For example, in 
the ISM, the turbulence is generated by supernova 
remnants and stellar-wind collisions, so one might 
expect the longest wavelength to be on the order of 
several parsecs or less.) Also, the lower end proba- 
bly corresponds to the scale at which the magnetic 
field couples most effectively to the particles, i.e., 
on the order of several gyration radii, since this is 
where the magnetic field loses most of its energy. 

Given the complexities and uncertainties in the 
global properties of the magnetic field structure 
within GMCs, we make several simplifying as- 
sumptions throughout this baseline work. Specif- 
ically, we assume a homogeneous medium and 
that all MHD fluctuations propagate with a uni- 
form (Alfvenic) speed va = 1 km s _1 . Although 
much of our analysis is dimensionless and there- 
fore easily scaled, we adopt fiducial values when 
dimensionalizing our results. Specifically, we as- 
sume that magnetic fluctuations have a maximum 
wavelength of A max — 1 pc (essentially the typical 
distance between stellar wind sources, as noted 
above). Further, we consider both the case of a 
purely turbulent field and the case of an underly- 
ing uniform magnetic field with a strong turbulent 
component. For the former case, we assume that 
the energy density of the turbulent field is equal 
to that of a 10 fiG uniform field. For the latter, 
we assume that the underlying uniform field has 
a magnetic strength of Bq = 10 /iG, and that the 
turbulent component has the same energy density 
as the uniform field. 

3. The Turbulent Magnetic Field 

A novel numerical method for analyzing the 
fundamental physics of ionic motion in a static tur- 
bulent magnetic field was presented by Giacalone 
& Jokipii (1994), who showed that ions in com- 
plete 3D situations readily cross the resulting mag- 
netic field. We generalize this pioneering work 
by considering time-dependent fluctuations that 
propagate with a uniform speed va (as first at- 
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tempted in a different context by Fraschetti & 
Meiia 2008). Within this framework, the mag- 
netic field through which cosmic rays of mass m 
and charge q propagate is expressed in terms of 
the gyration frequency via the parameter J~2(r, t) = 
qB(r, t)/mc. The total field is then written as the 
sum of a static background component l~2b(r) and 
a fluctuating, time dependent component <5J7(r, t), 
but we note that it is not necessary to have a back- 
ground component, and for cases where such a 
component exists, fluctuations need not be small. 
Further, a time-dependent turbulent electric field 
5E(r, t) must also be present (as required by Fara- 
day's law; Fraschetti & Melia 2008). As shown be- 
low, SE << SB for molecular cloud environments 
and, as such, the effects of such an electric field 
may be ignored in the analysis presented here. 

The turbulent magnetic field is generated by 
summing over a large number N of randomly 
polarized transverse waves of wavelength A„ = 
27r/£;„: 

jv 

(517 (r, t) = fi„ [cos a n y' ± i sin a n z'] 

n=l 

exp [ik n {x' - v A t) + i/3„] , (7) 

where k x = k mm = 2n/X max and k N = k max = 
27r/A m i n are, respectively, the wavenumbers cor- 
responding to the maximum and minimum wave- 
lengths associated with the turbulent field, the an- 
gle a n and phase (3 n are randomly selected be- 
tween and 27r, and the random choice of ± se- 
lects the helicity of the wavevector about the x' 
axis. The corresponding turbulent electric field is 
given by 



N 



<5E(r, t) — — fl n [±i sin a n y' — cos a n z'] 

q c f— ; 

n— 1 

exp [ikn(x' - VAt) + ij3 n ] ■ ( 8 ) 

The determination of the random polarization 
of each wavevector k n in the laboratory frame is 
accomplished via the two-angle rotation matrix 



R = 



cos 9 n — sin 9 n cos 4> n 
cos 6 n cos 4>n 
sin 4> n 




sin 8 n sin <f) n 
— cos 9 n sin 4> n 

COS 4>n 



0) 



where < tfi n < 2ir, and < cos^ n < 1 
are selected randomly (for a total of five ran- 



dom components for each value of n) Through- 
out this work, the turbulent field structure at 
any position r is calculated by summing over 
N = 25 log 10 [A max /A m in] values of wavevectors 
k n , evenly spaced on a logarithmic scale between 
k m i n and k max (as justified in §5). Specifically, the 
particle position in the primed frame r' = R • r is 
used to calculate the real part of the turbulent 
magnetic field for each wavevector k n , as given by 



Re{5n(T,t)' n } = 
tt n {cosa n cos [k„. (x' - v A t) + j3 n ] y 
±sina n sin [k n {x 1 - v A t) + /3„] z'} . 



(10) 



Since each k n component is randomly oriented 
(i.e., has its unique value of y' and z'), one must 
perform the rotation back to the unprimed frame 
Sn(r) k = R • Sn(r)' k (where R R = I— e.g., 
= Rj t i) before performing the sum over n. 
The desired spectrum of the turbulent magnetic 
field is set through the appropriate choice of T in 
the scaling 



SI 



k n 



-r 



Afc„ 
Afci 



-r+i 



(11) 



(as we have indicated, T = 1 for Bohm, 3/2 for 
Kraichnan, and 5 /3 Kolmogorov) , where the quan- 
tity Jli is set by a parameter £ that specifies the 
energy density of the turbulent field via the defi- 
nition 

"Jfe l~ r+1 

g =W 2 . (12) 



We note that for our adopted scheme, the value of 
Ak n /k n is the same for all values of n. We further 
note that £ = 2 corresponds to the real part of the 
turbulent field having the same energy density as 
a uniform field fi since Sft ■ Stt* = 2Re{Sft} 2 . 
Here we assume that there are a sufficiently large 
number of randomly polarized transverse waves so 
that the cross terms of the above dot product can- 
cel each other out. 

4. Dimensionless Equations of Motion 

The equations that govern the motion of rel- 
ativistic charged particles through the turbulent 



1 The ZX rotation scheme adopted here differs from that 
presented in Giacalone & Jokipii (1994). 
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medium are 



and 



dt mc \ 7 



dr 



dt 



(13) 



(14) 



where u = 7v/c and 7 is the particle Lorentz fac- 
tor. As can be seen from the form of Equations (7) 
and (8), SE~(v A /c)SB. Since MHD fluctuations 
in molecular clouds are expected to propagate at 
speeds of va ~ 1 km s -1 , SE « SB, and the 
electric field has a negligible effect on the local 
particle motion for particle speeds approaching c. 
However, electric fluctuations can significantly ac- 
celerate charged particles given a sufficiently long 
time (Fraschetti & Melia 2008). Under the most 
ideal conditions, turbulent fields can energize pro- 
tons in a time At by an amount 



AE„ = e8E cAt « e SB v A At . 



(15) 



Such an ideal acceleration, however, can only oc- 
cur for time intervals At < A max /c. For the pa- 
rameter values adopted here (v A = 1 km s _1 , 
SB = 10 yuG, A max = 1 pc), this ideal acceleration 
may only last for ^3 yrs and energize particles 
by an amount AE p w 0.03 TeV. For longer time 
intervals, the process becomes stochastic and the 
particle energy increases as AE p oc \ft. A reason- 
able upper limit to the increase in particle energy 
as a function of time is therefore given by 



AE. 



p; max 



0.01 TeV 



lyr) 



1/2 



(16) 



In order to both confirm this result and to ob- 



tain a more exact value for AE V 



we have 



J p\ max; 

solved Equations (13) and (14) for protons moving 
trough a turbulent field characterized by T = 3/2, 
Amax = 1 pc, A min = 10~ 4 pc, an energy density 
equal to that of a uniform B = 10 pG magnetic 
field, and our adopted fiducial value of Va = 1 km 
s _1 . Since the focus of our paper is on relativistic 
particles whose radius of gyration 



Rg - 



~qBo~ 



1.08 x 10~ 4 pc 



E p 
ITeV 



Bp 
10 mG 



(17) 



falls within the values of A m ; n and A max , we have 
solved the resulting equations of motion for both 
a 10 2 TeV and a 10 3 TeV proton. The resulting 
change in energy |AZ? p | as a function of time for 
both particles is shown in Figure 1, and clearly 
demonstrates a random-walk behavior (for which 
|A_E p | oc \/t) with fluctuations superimposed. In 
addition, we find that Equation (16) — as repre- 
sented by the dashed line in Figure 1 — provides a 
good upper limit for |A_E p |. Since we focus our 
discussion on particle energies in excess of 1 TeV 
and diffusion times less than 10 4 years, this test 
calculation shows that we may justifiably ignore 
the effects of the electric field in our work. 

To simplify the analysis, we define a dimension- 
less time t = t/to, where t is the inverse of the 
gyration frequency multiplied by the Lorentz fac- 
tor for a particle with charge q = Ze and mass 
m = ArriH in a reference field B , as given by the 
expression 



t - 7 



3.5 x 10~ 4 yrs 



E P \ (A 
ITeV / \Z 



Bo 



10 ^G 



.(18) 



We also define a corresponding dimensionless ra- 
dius vector r = r/R g . Since we ignore the electric 
field SE, |u| = 7«/c is a constant of the motion. 
Thus, for relativistic particles (v w c), setting the 
value of R g also sets the value of to (and vice versa) 
since R g = cto. 

Ignoring the electric field, the equations of mo- 
tion for highly relativistic particles can then be 
written in dimensionless form as 



du „ - 
^ =UXB ' 



and 



dr 

d7 



u 



(19) 
(20) 



where B = B/jB and u = u/|u| 



5. The Case of a Purely Turbulent Field 

In our formalism, the trajectory of a particle 
moving through a purely turbulent field is fully 
described by the four dimensionless parameters T, 
ua = va/c, A min = A min /i? g and A max = A max /i? g 
(related to the rigidity p through the expression 
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Fig. 1. — The magnitude of the change in par- 
ticle energy |A£? P | as a function of time for pro- 
tons with initial energies of 10 2 TeV (solid curve) 
and 10 3 TeV (dotted curve) moving through tur- 
bulent magnetic and electric fields characterized 
by T = 3/2, A max = 1 pc, A min = 1(T 4 pc, and 
va = 1 km s _1 . The turbulent magnetic field has 
an energy density equal to that of a uniform 10 
/iG field. The dashed line represents the value of 
the upper limit |A.E p; max | given by the expression 
in equation (16). 



Amax = 27r/p), along with the adopted prescrip- 
tion for setting the N values of wavevectors k n 
discussed below Equation (9). It is important to 
note that as the particle moves through the field, 
the radius of gyration changes depending on the 
field strength being sampled. Within this context, 
Bq is taken to be the field strength of a uniform 
field whose energy density equals that of the tur- 
bulent field. In turn, the value of R g represents a 
characteristic value for a particle's radius of gyra- 
tion. 

We begin our analysis by considering how mo- 
tion through a time-dependent turbulent field dif- 
fers from that of a static turbulent field (va = 0). 
To this end, we calculate the trajectory of a par- 
ticle over a time T max = 10 5 A max for the case 
r = 3/2 and the following four sets of wavelength 
ranges [A min ,A max ]: [3,300]; [0.3,30]; [0.03,3]; and 
[0.003,0.3]. We plot the displacement f of each 
particle as a function of (the dimensionless) time 
r in Figure 2 for the case of a static field (va = 0), 
and in Figure 3 for the case of a timc-dcpcndcnt 
magnetic field with an adopted fiducial value va 
= 1 km s _1 . The long-dashed lines serve as a ref- 
erence and have slopes of 1/2. 

Figures 2 and 3 illustrate three important 
points. First, particles with a radius of gyration 
below the range of turbulent wavelengths may 
eventually get trapped in a static field, as can be 
seen by the fact that f is constant at times r > 10 5 
for the Xmin = 3 particle (solid line in Figure 2). 
To gain insight into this phenomenon, we plot in 
Figure 4 the dot product B • v as a function of 
time for the trapped particle shown in Figure 2 
(solid curve) . One sees that trapping occurs when 
particles move nearly perpendicular to the local 
magnetic field, oscillating in a sort of local mag- 
netic bottle. As can be seen in Figure 3 from the 
solid line at times r > 10 , time-dependent fluc- 
tuations will disrupt this trapping on an expected 
timescale tmhd ~ ^min/uA (~ 10 6 for the solid 
curve shown in Figure 3). Second, once parti- 
cles with radii of gyration smaller than A max have 
moved beyond a (dimensionless) distance ~ A max , 
their displacement scales as f oc r 1 / 2 . Finally, par- 
ticles with a radius of gyration greater than the 
maximum turbulent wavelength are not strongly 
affected by local turbulence. The motion of such 
(highly-energetic) particles will not be considered 
in our analysis. 
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Fig. 2. — The displacement f as a function of time 
r for four particles moving through a static (va = 
0) turbulent magnetic field with index F = 3/2. 
The values of A m ; n and A max correspond to the 
following curves: 3, 300 (solid); 0.3, 30 (short- 
dashed); 0.03, 3 (dotted); 0.003, 0.3 (dot-dashed). 
The long-dashed line serves as a reference and has 
a slope of 1/2. 
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Fig. 3. — The displacement f as a function of time 
t for four particles moving through a temporally 
fluctuating (va = 1 km s _1 ) turbulent magnetic 
field with index T = 3/2. The values of A m ; n 
and A max correspond to the following curves: 3, 
300 (solid); 0.3, 30 (short-dashed); 0.03, 3 (dot- 
ted); 0.003, 0.3 (dot-dashed). The long-dashed 
line serves as a reference and has a slope of 1/2. 



The motion of charged particles through a tur- 
bulent magnetic field is chaotic in nature. As such, 
a complete analysis requires a statistical approach. 
We have therefore performed a suite of experi- 
ments designed to adequately sample our parame- 
ter space. Specifically, each experiment is defined 
by a choice of the parameters T, A m i n , and A max . 
We adopt the value of va = 1 km s _1 , although in 
the absence of particle trapping, our results will 
not be sensitive to this chosen value. For each 
run, we calculate the trajectory of N p particles in- 
jected randomly from the origin for a time T maxi 
with each particle sampling its own unique mag- 
netic field structure (i.e., the values of a n , /?„, 9 n , 
(f> n and the choice of a ± are chosen randomly 
for each particle). The suite of experiments per- 
formed for the case of a purely turbulent field are 
summarized in Table 1. 

We plot the distributions of x = x/R g and 
f = r/Rg at time r = 10 3 A max for experiment 

2 in Figures 5-6. (The corresponding distribu- 
tions for experiments 1 and 3 are qualitatively 
very similar.) Since the particles at this time 
have fully sampled the turbulent structure of the 
field, the distributions of their positions x, y and 
z are expected to be normal. For a purely tur- 
bulent field, all three distributions are expected 
to have mean values of zero and equal variances 
(within the expected statistical fluctuations) . Fur- 
thermore, since motion along any axis is indepen- 
dent of the others, then the displacement vector 
f = y/x 2 + y 2 + z 2 has three independent orthog- 
onal components, each of which follow a standard 
normal distribution. As such, the f values should 
be distributed according to a chi distribution with 

3 degrees of freedom. To illustrate these points, 
we include the corresponding Gaussian curve de- 
rived from the mean and variance in Figure 5, and 
the corresponding k — 3 chi distribution in Fig- 
ure 6. As illustrated by our results, cosmic-ray 
diffusion through turbulent magnetic fields is well 
represented by Gaussian statistics. 

The median, mean and rms values of the f dis- 
tribution shown in Figure 6 are denoted, respec- 
tively, by the vertical dotted, short-dashed, and 
long-dashed lines. Although each of these out- 
put measures characterize the distribution, we will 
adopt the mean value (r) of the particle displace- 
ments as our primary output measure, and calcu- 
late its value at several times r for each experi- 
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Fig. 4. — The dot product between the field direc- 
tion and particle direction of motion as a function 
of t for the trapped particle shown in Figure 2 
(solid curve). 
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Fig. 5. — The distribution of x values at time 
t = 10 3 A max for experiment 2 (histogram), super- 
imposed with a Gaussian function (black curve) 
with the same mean and variance. 



ment performed (as listed in Table 1). In order 
to determine how sensitive the value of our out- 
put measure is on A m i n , we compare the results 
of experiments 1-3 with those of experiments 7, 
9 and 11 in Figure 7. As clearly illustrated by 
the overlap between the results from experiments 
1 (open triangle) and 7 (solid triangle), 2 (open 
square) and 9 (solid square), and 3 (open circle) 
and 11 (closed circle), particle diffusion depends 
primarily on the maximum turbulence wavelength 
Amax, and is not sensitive to the minimum tur- 
bulence wavelength A m ; n , so long as the radius of 
gyration is greater than the minimum turbulence 
wavelength (see discussion in §7). Our analysis 
is therefore greatly simplified in that there is only 
one primary parameter - A max - that dictates how 
particles diffuse through a purely turbulent field 
with a specified value of T. We also note that 
the values of (f) clearly exhibit the t 1 / 2 depen- 
dence associated with a diffusion process (though 
particles with R g <~ A max have motions interme- 
diary to their counterparts with smaller radii of 
gyration and the free-streaming motion of their 
counterparts with greater radii of gyration). 

We focus the rest of our analysis on cases for 
which the particle gyration radius falls comfort- 
ably within the range of the maximum and mini- 
mum turbulence wavelengths so that particles un- 
dergo actual diffusion — that is, for which A max 3> 
1 3> A m i n . To do so, we consider a turbulent field 
with a dynamic range in wavelengths that span 
either four or five orders of magnitude. We note, 
however, that the minimal dependence that par- 
ticle diffusion has on the smallest wavelength im- 
plies that our results can be extrapolated to lower 
values of A m ; n (see discussion in §7). 

A fundamental issue in this analysis is what 
value of N will allow our discrete treatment of the 
turbulent field to adequately represent a contin- 
uous field. Toward that end, we first note that 
the variance of the mean values of (f) is given by 
Vmean = crf/WN p . Based on the results presented 
in Figure 6, o> <~ (f)/2, so that the calculated 
mean of our sample population with N p = 200 is 
expected to be within 3<r mean — 1.5(f) / tJN p rts 
0.1(f) of the true (parent) value with <~ 99% con- 
fidence. We next perform experiments 6, 16 and 
25 with values of N = 50, 100, 200 and 300. The 
resulting values of (f) at time T max as a function of 
N are shown in Figure 8, where the error bars rep- 
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Fig. 6. — The distribution of f values at time 
t = 10 3 A max for experiment 2 (histogram), su- 
perimposed with a chi function of degree 3 (black 
curve) and scaled using the mean of the x, y and z 
distribution variances. The vertical dotted, short- 
dashed, and long-dashed lines represent the me- 
dian, mean and rms values for the distributions, 
respectively. 
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Fig. 7. — The value of (f) as a function of t for 
experiments 1 (open triangle), 2 (open square), 3 
(open circle), 7 (solid triangle), 9 (solid square) 
and 11 (solid circle). The dashed lines serve as a 
reference and have slopes of 1/2 and 1. 



resent the expected 3tr statistical error of 0.1(r). 
These results appear to justify our adoption of 
N = 25 log 10 [A max /A min ] presented in §2. 

The results of experiments 4-10, 12-18 and 19- 
25 are presented in Figures 9, 10 and 11, respec- 
tively. A self-similar pattern is clearly visible in 
these figures for cases with A max > 30, with a 
break in the slope of the curves from <~ 1 to 1/2 
occurring around t ~ A max /10 for T = 3/2 and 
5/3, and at r ~ 10 for T = 1. We note, how- 
ever, that the break is not smooth for the solid 
circles show in Figures 9, 11 and 12. These ir- 
regularities occur as particles with small radii or 
gyration make transitions from weakly perturbed 
propagation (for which (f) oc r) to diffusion (for 
which (f) oc T 1//2 ). This feature indicates that as 
particles with small radii of gyration make this 
transition after traveling a distance <~ 0.1A max , 
they are effectively "scattered" randomly in all di- 
rections, so that on average, their distance from 
the origin does not change appreciably until they 
truly reach the diffusion regime (i.e. they have 
been "scattered" numerous times). 

In order to put our results into a physical 
context, we consider relativistic protons moving 
through a purely turbulent magnetic field for 
which A max = 1 pc, and dimensionalize the re- 
sults of experiments 4-10 accordingly through 



a proper choice of R g 



/v max /A 



We note 



that setting a common value of A max for ex- 
periments 4-25 also sets a common value of 
tmax = T max t = 10 2 A max /c. The results arc 
presented in Figure 12. As previously noted, the 
solutions are nearly self-similar for particles whose 
radius of gyration is R g < 0.03A max . 

To better understand how a particle's gyration 
radius helps determine the nature of its motion, 
we plot in Figure 13 particle trajectories of three 
particles with different radii of gyration, each in- 
jected with identical velocity from the origin into 
the same turbulent (but static) magnetic field de- 
fined by T = 3/2, A max = 1 pc, A min = 10~ 4 
pc, and Bo = 10 ^G. The field line that passes 
through the origin is depicted by the thin black 
line. Particle trajectories are depicted by the blue 
(R g = 0.001 pc), green (R g = 0.01 pc) and red 
(R g = 0.1 pc) curves. Clearly, the nature of parti- 
cle motion differs for particles with R g < 0.01A max 
and R g > 0.01A max . For the former, particles are 
strongly coupled to field lines and their motion is 
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Fig. 8. — The value of (f) at T max as a function of 
TV for experiments 6, 16, and 25. The error bars 
represent the expected 3cr statistical error of 10%. 



Fig. 10. — The value of (f) as a function of r for 
experiments 12-18, for which T = 1. The dashed 
lines serve as a reference and have slopes of 1/2 
and 1. 
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Fig. 9. — The value of (f) as a function of r for 
experiments 4-10, for which T = 3/2. The dashed 
lines serve as a reference and have slopes of 1/2 
and 1. 
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Fig. 11. — The value of (f) as a function of r 
for experiments 19-25, for which T = 5/3. The 
dashed lines serve as a reference and have slopes 
of 1/2 and 1. 
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directly tied to the field line structure, whereas 
for the latter, particles "random walk" through 
the field. That is not to say that particles with 
small radii of gyration move smoothly along field 
lines. Rather, although they are scattered by the 
turbulent magnetic fields according to their ener- 
gies, their spread due to scatter is small compared 
to how far they propagate in the direction of the 
field. 

A central aspect of this work is a determina- 
tion of the relation between particle diffusion and 
energy. To that end, we define a dimensionless 
energy e = E p /Eq, where 



Eq — A max e Bq — 
9.2 x 10 3 TeV 



Bp 

lpc J V 10//G 



(21) 



which then yields the relation e = Zi?, 5 /A max = 
z KLt- We P lot tne values of (r)/A max at r max 
as a function of e in Figure 14 for experiments 
4-10, 12-18, and 19-25. Each set of results for 
a given value of T demonstrates a clear break at 
tb ~ 0.005, corresponding to particles with gy- 
ration radii R g ~ 0.005A max /.Z'. There is clearly 
a stronger dependence between (7-) and e above 
the break, presumably due to the fact that par- 
ticles with e « e b are strongly coupled to the 
field lines, as shown in Figure 13. As such, their 
diffusion is dictated primarily by the field struc- 
ture, and hence, becomes less sensitive to their 
energy/radius of gyration. Specifically, particles 
with radii of gyration smaller than ~ 0.005A ma a: 
will effectively scatter off field fluctuations that 
have a similar length scale as their gyration radius. 
In contrast, particles with sufficiently large gyra- 
tion radii effectively decouple from the field-lines 
(as is illustrated in Figure 13), and essentially ran- 
dom walk through the field on length scales equal 
to their gyration radius. Their motion, therefore, 
is not very sensitive to the nature of the small-scale 
fluctuations, as can be seen by the convergence of 
the output values in this regime for the r = 5/3 
and r = 1/2 cases. 

In order to put our results into a useful for- 
mat, we note that in the standard theory for par- 
ticle diffusion, the turbulent field index T is related 
to the diffusion coefficient index S (as defined in 
Equation 2) through the expression 6 = 2 — T. As 
such, the diffusion radius R diS oc E\~ V I 2 t 1 ! 2 . In 
turn, we express the particle diffusion length as a 
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Fig. 12. — The value of (r) as a function of time 
for experiments 4-10, for which V = 3/2. The 
results of these experiments are dimensionalized 
by assuming that A max = 1 pc for each case, and 
setting the value of R g accordingly. The dashed 
lines serve as a reference and have slopes of 1/2 
and 1. 




Fig. 13. — Trajectories of three particles injected 
with identical velocities from the origin into the 
same turbulent (but static) magnetic field, defined 
by T = 3/2, A max = 1 pc, and A min = 10~ 4 
pc. The colored curves denote the path of par- 
ticles with gyration radii 0.001 pc (blue), 0.01 pc 
(green), and 0.1 pc (red). The black curve denotes 
the magnetic field line passing through the origin. 
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function of energy and time through the expres- 



sion 



where 



(r) = K 



1/2 



tc 



3.3 yrs 



1 pc 



(22) 



(23) 



We then fit the three lowest-energy data points for 
each case shown in Figure 14 at time t — 100 t c , 
as illustrated by the dashed (r = 1), solid (r = 
3/2) and dash-dotted (r = 5/3) lines, where the 
corresponding values of A and a are given in Table 
3 for each value of T. In all cases, good fits are 
obtained with a = 1 - T/2 for e < 0.005. 

6. Uniform Field with a Turbulent Com- 
ponent 

We next consider a molecular cloud environ- 
ment threaded by a uniform magnetic field with 
a strong turbulent component. Specifically, we 
assume a magnetic field of the form B(r, t) = 
Bqz + <5B(r,i). In our formalism, the motion of 
a particle moving through such a field is then de- 
scribed by five dimensionless parameters: T, ua, 
Amax and rj. 

Observations of molecular clouds suggest that 
the magnetic fluctuations have amplitudes 5B ~ 
£?o. This finding follows from considering the 
observed non-thermal line-widths in molecular 
clouds (Larson 1981; Myers et al 1991) to result 
from MHD waves (e.g., Fatuzzo & Adams 1993; 
McKcc & Zweibel 1995; see Fatuzzo & Adams 
2002 for further discussion). We therefore con- 
sider the case that the magnetic energy density of 
the turbulent field equals that of the underlying 
field, thereby setting £ = 2 for all cases explored. 
The suite of experiments performed are summa- 
rized in Table 2. 

The introduction of the field B z has broken 
the isotropy, so we now plot both the distribu- 
tion of x = x/Rg and that of z = z/R g at time 
r = 10 2 A max for experiment 5 (Table 2) in Fig- 
ures 15 and 16, were the solid curves depict the 
corresponding Gaussians derived from the mean 
and variance of each distribution. As illustrated 
by our results, cosmic-ray diffusion through uni- 
form magnetic fields with strong turbulent compo- 



nents is fairly well represented by Gaussian statis- 
tics. 

The rms values of the particle positions x and z 
at several times t for each experiment are shown 
in Figures 17 and 18. As found for the purely 
turbulent field discussed in §5, the curves appear 
to be nearly self-similar, with a break in the slope 
of the curves occurring at around t ~ A max /10. 
Not surprisingly, particles diffuse further along the 
direction of the uniform field than they do across 
the field, with z rms ~ 5x rms . 

As noted in §5, a central aspect of this work is a 
determination of the relation between particle dif- 
fusion and energy. To that end, we plot the values 
of x rms / X max and z rms /X max at T max as a function 
of e in Figure 19 for experiments 1-21 listed in Ta- 
ble 2. In all cases, the data for diffusion along the 
underlying magnetic field direction is well-fit by a 
line. Likewise, the data for the diffusion across the 
underlying magnetic field is well-fit by a line for 
T = 1 and T = 3/2, but does exhibit a break at 
e ~ 0.01 for T = 5/3. 

Following the analysis presented in §5, we ex- 
press the particle diffusion lengths across and 
along the underlying uniform magnetic field 
through the expressions 



and 



= A n 



-A, 



E~ 



1/2 



1/2 



(24) 



(25) 



We fit the data in Figure 19 at time t = 100 t c , 
as illustrated by the dashed (r = 1), solid (r = 
3/2) and dash-dotted (r = 5/3) lines, where the 
corresponding values of A and a are given in Table 
4 for each value of T. In all cases except for x rms 
when r = 1, good fits are obtained with a = 1 — 
T/2 for the entire range of e explored. 

7. Comparison to Previous Work 

The transport properties for charged particles 
moving through turbulent magnetic fields was an- 
alyzed by Casse et al. (2002) using a method 
similar to that adopted in our work. Specifically, 
these authors performed extensive numerical ex- 
periments using the formalism developed by Gi- 
acalone & Jokipii (1994) in order to determine 
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Fig. 14. — The value of (r)/A max evaluated at 
Tmax as a function of e for experiments 4-10 
(r = 3/2), 12-18 (r = 1), and 19-25 (T = 5/3). 
The dot-dashed (T = 5/3), solid (r = 3/2) and 
dashed (T = 1) curves represent power-law fits to 
the data, as discussed in the text. 
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Fig. 16. — The distribution of z values at time 
r = 10 2 A max for experiment 5 in Table 2 (his- 
togram), superimposed with a Gaussian function 
(black curve) with the same mean and variance. 




Fig. 15. — The distribution of x values at time 
r = 10 2 A max for experiment 5 in Table 2 (his- 
togram), superimposed with a Gaussian function 
(black curve) with the same mean and variance. 



Fig. 17. — The value of x rms as a function of r 
for experiments 1-7 listed in Table 2. The dashed 
lines serve as a reference and have slopes of 1/2 
and 1. 
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^rms/Amax (open) evaluated at r maa; as a function 
of e for experiments 1-7 (r = 3/2), 8-14 (r = 1), 
and 15 21 (r = 5/3) listed in Table 2. The dot- 
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cussed in the text. 



the pitch angle, scattering rate, and the parallel 
and perpendicular spatial diffusion coefficients for 
a wide range of rigidities and turbulence levels. 
Both parallel and perpendicular diffusion coeffi- 
cients are plotted versus rigidity p = R g k min = 
2TT/X rn ax for several different values of turbulence 
level 



(SB' 



V = 



(26) 



We note that Casse et al. (2002) employed 
two different methods to construct their magnetic 
fields. For 77 = 1 (which represents a purely turbu- 
lent field), these authors adopted the same scheme 
presented in our work, and used a dynamic range 
in wavelengths of A max /A min = 10 4 . For all other 
cases, the magnetic field was constructed using a 
fast-Fourier transform (FFT) algorithm to set up 
the magnetic field on a discrete grid in configu- 
ration space. An interpolation scheme was then 
used to calculate the field at any point in space. 
For this latter method, X max /\ m i n = 128 for most 
cases. 

Our analysis extends the work of Casse et al. 
(2002) in two ways. First, while these authors 
focused exclusively on Kolmogorov diffusion, we 
also consider Bohm and Kraichnan diffusion. Sec- 
ond, we extend considerably the dynamic range 
of turbulence wavelengths, especially for the case 
of a uniform field with underlying turbulence. In 
addition, we focus our results to the propagation 
of cosmic-rays in molecular cloud environments. 
Nevertheless, sufficient overlap exists for a direct 
comparison of a subset of our works. Specifically, 
experiments 19 - 25 listed in Table 1 (purely turbu- 
lent field) can be compared directly with the 77 = 1 
data presented in Figure 4 of Casse et al. (2002). 
In order to do so, we calculate the corresponding 
diffusion coefficients 



D 

RgC 



(A3: 2 ) 



(27) 



where Ax is the particle displacement (from the 
origin) along the x direction (although all direc- 
tions are equivalent) evaluated at time r = T max . 
We note that this method, while not exactly sim- 
ilar, is analogous to that adopted by Casse et al. 
(2002). As shown in Figure 20, our results (open 
squares) are in agreement with those of our pre- 
decessors (filled circles), with the dotted line de- 
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noting the value of p m in — 27r/A max used in their 
calculations. 

We next compare our results from §6 for the 
case of a uniform field with underlying turbulence 
to the r\ = 0.46 case presented in Figures 4 and 5 
of Casse et al. (2002). To do so, we calculate both 
perpendicular and parallel diffusion coefficients 



D ± _ (Ax 2 ) 



RgC 



2t„ 



RgC 



(Az 2 ) 



(28) 



for experiments 15 - 21 in Table 2. We com- 
pare our results (open squares and circles) to those 
of our predecessors (filled squares and circles) in 
Figure 21. The dotted line denotes the value of 
Pmin = 2Tr/\ max used by Casse et al. (2002) for 
this case. As expected, the results are in good 
agreement for p > p m i ni but deviate for lower val- 
ues of rigidity, further illustrating our conclusion 
from §5 that particle diffusion is not sensitive to 
the value of Xmin so long as R g > X m in- 

8. Application to Cosmic-ray Diffusion in 
Molecular Clouds 

One of the original motivations for this calcula- 
tion was to determine what kind of injection pro- 
file would be required in order to correctly inter- 
pret the apparent correlation between the diffuse 
7-ray emissivity and the distribution of molecular 
gas in the interstellar medium. Such a correlation 
between 7-ray intensity maps and the large-scale 
features of the diffuse gas was first noted in ob- 
servations (E-y > 100 MeV) with the SAS-2 and 
COS B satellite telescopes, combined with radio 
data that reveal the column density of interstellar 
hydrogen. Later observations associated at least 
ten EGRET sources with SNRs expanding into 
MCs (Esposito et al. 1996; Combi et al. 1998, 
2001; Torres et al. 2003). More recently— and 
more spectacularly — a strong correlation between 
TeV emission and the molecular gas distribution 
at the Galactic center was demonstrated by HESS 
(Aharonian et al. 2006a; Wommer et al. 2008). 
These data lend support to the idea that the low 
latitude 7-ray emission is mainly due to the de- 
cay of neutral pions produced by the scattering of 
cosmic rays with protons in the ambient medium 
rather than from bremsstrahlung or inverse Comp- 
ton (IC) scattering. 

In their assessment of this effect, Aharonian 
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Fig. 20. — Comparison of the diffusion coefficients 
calculated for experiments 19-25 listed in Table 
1 (open squares) with the corresponding diffusion 
coefficients presented in Figure 4 of Casse et al. 
2002 (filled circles). The dotted vertical line de- 
notes the value of p m in = ^/K nax adopted by 
this earlier work for the results shown here. 
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sion coefficients calculated for experiments 15 - 
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Pmin = 2Tr/X max adopted by this earlier work for 
the results shown here. 
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and Atoyan (1996) argued that the principal re- 
gion of interest for the 7r°-decay 7-ray emission 
ought to lie within an R < 100 pc region surround- 
ing the cosmic-ray source. Within this distance 
of a "typical" particle accelerator, a total energy 
output of W p <~ 10 50 erg translates into a mean 
particle energy density of w p = W p / (4/3)irR 3 s=y 
0.55(VK P /10 50 erg)(i?/100 pc)" 3 eV/cm 3 , which 
may significantly exceed the average level of the 
"sea" of galactic cosmic rays with energy den- 
sity wo ~ 1 eV/cm 3 . Therefore, in a 1°-10° re- 
gion around a cosmic-ray source (depending on 
the distance to the source), we should expect to see 
higher than average 7-ray emission. In addition, if 
the diffusive propagation of cosmic rays is energy- 
dependent, the resulting 7-ray spectrum will dif- 
fer from the 7-ray spectrum produced by galactic 
cosmic rays (e.g., Fujita et al. 2009). Thus, the 
possibility of having several dense giant molecu- 
lar clouds (GMCs) in close proximity to a parti- 
cle accelerator will not only produce higher than 
average levels of 7-rays but may give the appear- 
ance that there are multiple distinct cosmic-ray 
sources or, due to the limited angular resolution of 
instruments like EGRET, an extended cosmic-ray 
source. Accurately predicting the spatial and tem- 
poral evolution of the 7-ray spectrum produced 
by a particle accelerator may therefore lead to 
the classification of tens of unidentified EGRET 
sources. 

In order to apply our results from §§5 and 6 
to molecular cloud environments, we consider the 
ideal case of a single impulsive cosmic-ray source 
surrounded by a homogeneous molecular cloud of 
radius R. While the value of A max is not known 
for such environments, one would expect its value 
to be constrained from below by the size of dense 
cores (<~ 0.1 pc) and from above by the size of 
the actual cloud (~ 10 — 20 pc). We there- 
fore adopt the intermediary value of A max = 1 
pc in our discussion (although we keep \ ma x in 
our scaled equations below). The energy range 
10~ 4 < e < 0.1 of our work (as shown in Figures 
14 and 19) thus corresponds to a true particle en- 
ergy range of 1 < E p < 10 3 TeV. In turn, since 
only ~ 10% of a relativistic protons' energy goes 
into the tt photon decay channel for pp scattering 
(see, e.g., Fatuzzo et al. 2006), the corresponding 
energy range of 7-rays resulting from the inter- 
action of these CR's and the ambient molecular 



cloud medium is 0.1 < e 7 < 10 2 TeV, which falls 
within the range observable by HESS. 

As shown by Aharonian and Atoyan (1996), the 
energy loss rate of protons with energies needed to 
produce 7r°-decay 7-rays is dominated by nuclear 
energy losses due to pp scattering with the ambi- 
ent medium. The lifetime of the protons, r pp , de- 
pends on the pp-scattering cross-section, a ppi and 
the inelasticity parameter, k. Over a broad range 
of proton energies, neither of these quantities sig- 
nificantly varies so the usual method is to adopt 
the constant average values a pp w 40 mb and n w 
0.45 (see, e.g., Markoff et al. 1997). That being 
the case, the proton lifetime becomes independent 
of proton energy: 



t pp = (ncK(Tp P ) 1 « 3 x 10 5 yr ^ 



n H2 



100 cm" 3 



(29) 

where n is the number density of ambient protons 
(i.e., n = 2n H2 )- 

We compare this timescale to the particle es- 
cape time r e , defined here as the time it takes CR's 
to diffuse a distance (r) = R for purely turbulent 
fields, and (z rms ) = R if an underlying uniform 
field threads the molecular cloud. For the inter- 
mediary case of Kraichnan diffusion, Equations 21 
- 23 can be combined to yield the expression 



T e;t urb ~ 4 x 10 5 yrs 



R 



20 pc 



En 



ITeV 



-1/2 



1 pc 



-1/2 



Bp 
10 mG 



1/2 



• (30) 



Likewise, Equations 21, 23 and 25 can be com- 
bined to yield the expression 



En 



ITeV 



T e; unif ~ 10 yrs 



-1/2 / . N -1/2 



R 



20 pc 



1 pc 



Bp 
10 mG 



1/2 



• (31) 



As suggested by Figure 1, injected particles will 
therefore gain a modest energy of ~ 1 TeV due 
to acceleration from the turbulent electric fields 
before they escape. As such, the fits to the data 
presented in Figures 14 and 19 (as summarized 
by Equations 22 - 25) cannot be extrapolated to 
lower energies for molecular cloud environments 
(since e = 10~ 4 represents a particle energy of 
0.92 TeV under the assumed conditions). 
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Given the similarity between t pp and r e , a sig- 
nificant fraction of > TeV CR's will likely undergo 
pp scattering before escaping from the molecu- 
lar cloud environment. As this fraction decreases 
with increasing energy, the resulting 7-ray spec- 
trum will be softer than that of the injected par- 
ticle spectrum (e.g., Fujita et al. 2009). In ad- 
dition, if magnetic fields in molecular clouds are 
purely turbulent, then the break in the (r) - e 
data shown in Figure 14 at e ~ 0.005 - which for 
the assumed conditions corresponds to a value of 
E p ~ 50 TeV - would likely produce a break in an 
observed 7-ray spectrum at around e 7 ~ 5 TeV. 
Such a break would not be observed if molecu- 
lar clouds are threaded by an underlying uniform 
magnetic field (see, e.g., Figure 19). 

The total 7-ray luminosity expected from our 
assumed molecular cloud with a singe injection 
source W p is independent of the escape time, as 
can be seen through the simple estimate 



'PP/ \ 'e 



10 36 ergs" 
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(32) 



where / takes into account that only <~ 10% of 
the relativistic protons' energy goes into the tt 
photon decay channel. 

Interestingly, this value is in reasonable agree- 
ment with the w 10 35 erg s _1 luminosities in the 
0.1 - 100 GcV band inferred for four SNR's inter- 
acting with molecular clouds (G349.7+0.2; CTB 
37A; 3C 391; G8.7-0.1) observed by the Large 
Area Telescope on board the Fermi Gamma-ray 
Space Telescope (Castro & Slane 2010). Two of 
these SNRs (CTB 37A and G8.7-0.1) are also pos- 
sible counterparts to HESS sources with implied 
luminosities in the 0.2 - 10 TeV band of w 5 x 10 34 
ergs s _1 and k 2 x 10 35 erg s _1 , respectively, 
and three additional HESS sources coincident with 
SNRS G338.3-0.0, G12. 82-0.02 and W41 have im- 
plied 0.2 - 10 TeV luminosities of w 2 x 10 35 erg 
s _1 , 3 x 10 34 erg s _1 , and 4 x 10 34 erg s _1 , re- 
spectively (Aharonian et al. 2006b) . Finally, we 
note that the above expected luminosity also falls 
within the range 1 x 10 34 -4 x 10 36 ergs s _1 inferred 
from observations of the EGRET SNRs, although 
the energy range of this instrument only goes up 
to - 30 GeV. 



9. Conclusion 

We have investigated how high-energy CR's 
propagate through molecular cloud environments 
using a modified numerically based formalism de- 
veloped by Giacalone & Jokipii (1994) for the gen- 
eral study of cosmic-ray diffusion, thereby pro- 
viding a baseline analysis for two magnetic field 
configurations: 1) a purely turbulent field; and 
2) a uniform magnetic field with a strong turbu- 
lent component. We have focused most of our 
analysis on cases for which the particle gyration 
radius R g falls comfortably within the range of 
wavelengths shaping the turbulence. For a purely 
turbulent field, the trajectory of a particle is fully 
described by four dimensionless parameters. How- 
ever, we have found that the diffusion of an ensem- 
ble of particles through a turbulent field (char- 
acterized by the index V) depends primarily on 
only one of these — the dimensionless scale length 
Amax = A max /i? g . For a uniform field with a tur- 
bulent component, CR diffusion depends on one 
additional dimensionless parameter — the ratio of 
turbulent field energy density to the uniform field 
energy density. 

Given the chaotic nature of particle motion 
through turbulent magnetic fields, we performed 
a suite of statistical experiments as defined by 
the dimensionless parameters listed in Table 1 (a 
purely turbulent field) and Table 2 (a uniform field 
plus a strong turbulent component). Specifically, 
we calculated the trajectory of N p particles in- 
jected randomly from the origin for a time T max 
for each experiment, with each particle sampling 
its own unique (and randomly selected) magnetic 
field structure. The resulting distributions of par- 
ticle displacement along a given axis were found 
to be well described by Gaussian profiles with the 
same mean and variance, thereby justifying our 
use of the mean of the particle displacements (r) 
as our output measure for characterizing the dif- 
fusion of particles through purely turbulent fields, 
and the rms values of the particle positions x rms 
and z rms as our output measures for a uniform 
magnetic field Bqz with a strong turbulent com- 
ponent. We have found that after an initial time 
during which particles travel a distance ~ A max , 
each of these output measures scales as \ft, as ex- 
pected for a diffusion process. 

The results of our analysis indicate that particle 
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diffusion behaves differently for gyro radii in the 
ranges 0.01A max < R g < A max and R g < 0.01A max . 
Specifically, we have found that in the former 
case, particles "random walk" through the field, 
whereas for the latter, particles are strongly cou- 
pled to field lines and their motion is directly tied 
to the field line structure. In turn, the distance 
over how far particles diffuse in purely turbulent 
fields as a function of energy exhibits a clear break 
at the point where the particle's gyration radius 
R g « 0.005A max . 

Comparing our results with those obtained in 
earlier works, we find good agreement with pre- 
vious results obtained using the same formalism 
(e.g., Casse et al. 2002). In addition, our re- 
sults are well-fit by the "standard" scaling law 
l— r /2 

i?diff E p often invoked in the literature. We 
provide simple scaling relations between mean dif- 
fusion lengths and energy for both magnetic field 
profiles considered. We note, however, that these 
scaling-laws lead to a significant underestimation 
of the diffusion lengths for the case of purely tur- 
bulent fields at energies E p > 0.005 X max e Bq. In 
addition, the index 1 — T/2 is not valid for the 
case of Bohm diffusion (r = 1) perpendicular to 
an underlying uniform magnetic field. 

The results of our work have important con- 
sequences for properly connecting 7-ray spectra 
associated with molecular clouds to the underly- 
ing particle populations. We find that a signifi- 
cant fraction of > TeV CR's will likely undergo pp 
scattering before diffusing out of a molecular cloud 
environment. As this fraction decreases with in- 
creasing energy, the resulting 7-ray spectrum will 
be softer than that of the injected particle spec- 
trum (e.g., Fujita et al. 2009). In addition, if 
magnetic fields in molecular clouds are purely tur- 
bulent, then the break in the (r) - e dependence 
(as shown in Figure 14) is expected to produce 
a corresponding break in an observed 7-ray spec- 
trum at around e 7 ~ 5 TeV. Such a break would 
not be observed if molecular clouds are threaded 
by an underlying uniform magnetic field. 

The work we have reported here has conse- 
quences for other types of high-energy sources as 
well. For example, the compact object IE 1740, 
embedded within a molecular cloud at the galactic 
center, produces a jet of (presumably) relativistic 
electrons and positrons (Misra & Melia 1993) that 
eventually diffuse into the surrounding medium. 



The diffuse radio inensity from this region provides 
some measure of the lepton injection rate, but it 
clearly also depends on the energy-dependent dif- 
fusion rate through the molecular gas. The results 
reported here for proton diffusion cannot be di- 
rectly generalized to the case of positrons, but we 
anticipate seeing qualitative similarities between 
the two once we have completed the analogous 
positron simulations. 

The galactic center hosts a complex array of dif- 
fuse emitters, in addition to the TeV sources we 
have discussed in this paper. A proper analysis 
of the underlying nonthermal particle population 
producing this emission should therefore include 
observations at 7-ray (and even hard X-ray) ener- 
gies, in addition to the HESS data we have con- 
sidered here (see, e.g., Belanger et al. 2004; Rock- 
efeller et al. 2004). In future work, we will more 
closely examine the observational consequences of 
the different behavior of CR's above and below 
the break energy Ef,, particularly as it impacts the 
diffuse broadband emission within ~20 pes of the 
supermassive black hole Sgr A*. 

Of course, Sgr A* itself is apparently a signif- 
icant accelerator of relativistic electrons and pro- 
tons (Liu et al. 2006), the latter diffusing (Ballan- 
tyne et al. 2007) through the captured, accreting 
gas (Ruffert and Melia 1994; Falcke et al. 1997) 
into the surrounding medium, possibly producing 
the HESS point source coincident with the black 
hole. However, attempts at reconciling this TeV 
emission with the longer wavelength radiation pro- 
duced closer to the center have been hampered by 
the uncertain energy- dependence of this diffusion 
process. As we have discussed in this paper, a de- 
tailed knowledge of the diffusion coefficient is es- 
sential for meaningfully connecting the observed 
spectrum to the underlying nonthermal particle 
population. We will be applying the conclusions 
reached here to this important problem and will 
report the results elsewhere. 
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Table 1 

Experiments for A Purely Turbulent Field 
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Table 2 

Experiments for A Uniform Field Plus Strong Turbulence 





1 


V 


^min 


\ 

^max 


AT 

N p 


'max i ^max 


1 


3/2 


2 


0.1 


10,000 


200 


10 2 


2 


3/2 


2 


0.3 


3,000 


200 


10 2 


3 


3/2 


2 


0.1 


1,000 


200 


10 2 


4 


3/2 


2 


0.03 


300 


200 


10 2 


5 


3/2 


2 


0.01 


100 


200 


10 2 


6 


3/2 


2 


0.003 


30 


200 


10 2 


7 


3/2 


2 


0.001 


10 


200 


10 2 


8 


1 


2 


0.1 


10,000 


200 


10 2 


9 


1 


2 


0.3 


3,000 


200 


10 2 


10 


1 


2 


0.1 


1,000 


200 


10 2 


11 


1 


2 


0.03 


300 


200 


10 2 


12 


1 


2 


0.01 


100 


200 


10 2 


13 


1 


2 


0.003 


30 


200 


10 2 


14 


1 


2 


0.001 


10 


200 


10 2 


15 


5/3 


2 


0.1 


10,000 


200 


10 2 


16 


5/3 


2 


0.3 


3,000 


200 


10 2 


17 


5/3 


2 


0.1 


1,000 


200 


10 2 


18 


5/3 


2 


0.03 


300 


200 


10 2 


19 


5/3 


2 


0.01 


100 


200 


10 2 


20 


5/3 


2 


0.003 


30 


200 


10 2 


21 


5/3 


2 


0.001 


10 


200 


10 2 



Table 3 

Fitting parameters for Figure 14 



r 


A 


a 


l 
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0.56 
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Table 4 

Fitting parameters for Figure 19 
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